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1.0  INTRODUCTION 


The  effect  of  plasticity  on  the  rate  of  growth  of  fatigue  cracks  is 
significant  for  a  wide  range  of'  problems  associated  with  the  damage  tolerance 
assessment  of  aerospace  structures.  The  range  of  problems  includes  crack 
growth  from  cold  worked  fastener  holes,  crack  growth  through  plasticity  due  to 
local  notch  stresses,  crack  driving  force  for  thermal  gradient  fields  and 
welding  residual  strain  fields,  small  flaw  growth  in  high  nominal  stress 
fields,  and  numerous  related  problems.  These  problems  have,  of  course,  been 
analyzed  using  a  variety  of  approximate  analytical  or  numerical  procedures. 
However,  as  will  be  summarized  within  this  report,  many  of  these  earlier 
modeling  approaches  have  involved  errors  which  may  significantly  affect  the 
predicted  fatigue  crack  growth  life  of  the  structure.  The  current  research 
has  resulted  in  some  new  and  critical  insights  into  this  class  of  problems, 
while  providing  a  basis  for  improved  modeling  of  these  problems. 

The  current  research  makes  use  of  the  boundary  integral  equation  (BIE) 
method,  as  modified  to  account  exactly  for  the  elastic  crack  problem.  The 
usual  BIE  formulation  for  elastic  problems  reduces  the  numerical  problem  to 
one  of  modeling  the  boundary  data,  while  preserving  the  complete  interior 
solution  of  the  field  equations.  In  the  elastic  fracture  mechanics  problem, 
the  Green's  function  approach  is  used  wherein  the  BIE  is  modified  to  account 
for  the  presence  of  a  stress  free  crack  at  an  arbitrary  location  within  the 
structure.  The  use  of  the  Green's  function  for  the  crack  eliminates  the  need 


to  model  the  boundary  of  the  crack,  and  provides  a  complete  mathematical 
description  of  the  elastic  strain  field  within  the  body,  due  to  the  crack. 


This  clearly  contrasts  with  the  finite  element  method  which  requires  that  the 
crack  surface  and  the  interior  strains  be  modeled  with  some  set  of 


interpolation  functions. 

The  BIE  method  has  been  successfully  modified  to  account  for 
elastoplastic  response  by  a  number  of  investigators.  However,  extension  of 
the  fracture  mechanics  model  with  the  Green's  function  approach  has  not  been 
previously  demonstrated.  In  order  to  account  for  elastoplastic  response  with 
the  BIE  method  one  must  numerically  model  the  interior  plastic  strain  field. 

In  all  other  ways  the  elastoplasticity  solution  uses  the  standard  elastic  BIE 
formulation.  The  current  work  reports  on  the  successful  extension  of  the 
special  Green’s  function  formulation  for  the  fracture  mechanics  problem  to  the 
elastoplasticity  formulation.  Not  only  has  the  work  resulted  in  accurate 
models  of  crack  tip  plasticity  for  a  reference  problem,  but  it  has  shown  some 
important  new  analytical  and  numerical  results  for  cracks  growing  in  plastic 
strain  fields. 

The  second  year  of  the  contract  effort  focussed  on  the  crack  extension 
problem.  In  the  elastic  case,  a  direct  solution  method  for  fracture  mechanics 
weight  functions  was  established.  The  elastoplastic  problem  considered  the 
extension  of  the  elastic  crack  into  its  prior  plastic  wake.  The  effects  of 
crack  tip  overloads  on  retardation  or  acceleration  through  closure  and 
residual  stress  effects  are  included.  In  addition,  the  elastoplastic  BIE 
formulation  was  more  fully  exploited  for  problems  of  crack  growth  in  residual 
strain  fields  such  as  weldments. 

Some  work  addressed  improvements  in  a  new  flat  crack  BIE  formulation  for 
3D  fracture  mechanics  analysis.  The  majority  of  the  work  on  this  task  was 
funded  by  the  Internal  Research  Panel  of  Southwest  Research.  However,  some 


2.0  RESEARCH  OBJECTIVES 


Generally  speaking,  advanced  aerospace  structures  have  been  designed  for 
damage  tolerance  considerations  using  elastic  fracture  mechanics  models. 
Problems  associated  with  residual  plastic  strains  at  notches,  cold  worked 
fastener  holes,  weld  residual  strains,  and  thermal  gradient  loading  have  been 
modeled  using  elastic  superposition  methods  along  with  elastic  fracture 
mechanics  models.  Crack  tip  plasticity  is  involved  in  all  fatigue  crack 
growth  problems.  Crack  tip  plasticity  dominates  the  problem  of  predicting 
crack  growth  under  spectrum  loading  conditions  where  acceleration  and 
retardation  effects  are  important.  Finally,  the  small  flaw  problem,  wherein 
crack  growth  rate  is  apparently  accelerated  relative  to  the  large  flaw 
problem,  cannot  be  currently  explained  by  elastic  fracture  mechanics 
considerations. 

The  development  of  improved  models  for  the  crack  growth  problem  for  the 
full  range  of  these  problems  is  crucial  to  improved  damage  tolerance 
assessment  for  advanced  aerospace  structures.  The  overall  objective  for  the 
current  research  is  to  provide  a  new  basis  for  making  damage  tolerance 
assessments  through  numerical  modeling  of  crack  tip  behavior,  including  the 
effects  of  plastic  or  other  residual  strains.  The  elastoplastic  BIE  method  is 
the  basis  for  the  current  effort. 

*  The  first  goal  of  the  originally  proposed  program  was  to  extend  an 
existing  planar  elastic  fracture  mechanics  analysis  based  on  the  BIE 
methodology  to  the  analysis  of  plastic  zones  around  cracks.  The  second 
proposed  goal  was  to  establish  fundamental  results  for  crack  tip  elastoplastic 
behavior,  based  on  a  numerical  and  analytical  study  of  the  elastoplastic  BIE 
formulation.  The  third  proposed  goal  was  to  establish  the  credibility  of  the 


elastoplascic  BIE  formulation  relative  to  the  finite  element  method  for 
refined  numerical  analysis  of  the  nonlinear  fracture  mechanics  problem,  and  to 
apply  the  capability  to  important  problems  of  fatigue  crack  growth  modeling 
for  advanced  aerospace  structures.  The  goal  for  the  second  year  of  the  effort 
was  to  extend  the  research  to  the  problem  of  modeling  crack  extension  under 
elastoplastic  conditions. 

This  report  summarizes  key  findings  of  the  current  research  effort.  The 
next  section  summarizes  the  basic  two-dimensional  elastoplastic  formulation 
and  applications.  Included  in  this  work  are  the  preliminary  applications  of 
the  new  method  to  crack  extension  into  prior  plastic  zones.  The  next  section 
reports  on  the  use  of  the  new  3IE  formulation  for  elastic  crack  extension. 

This  new  result  allows  for  the  direct  computation  of  crack  weight  functions. 
The  last  section  reports  on  some  recent  work,  for  the  3D  BIE  fracture 
mechanics  formulation.  Some  contrast  with  the  2D  formulation  is  noted. 

Further  work  on  the  3D  problem  is  expected  in  the  subsequent  research  program. 


3.0  ELASTOPLASTIC  FRACTURE  MECHANICS  MODELING 


3 • 1  Review  of  the  Mathematics 

A  complete  treatment  of  the  elastic  formulation  for  the  Green's  function 
BIE  model  for  cracked  planar  problems  is  given  by  Snyder  and  Cruse  [ 1 ]  and 
Cruse  [2].  The  full  development  of  the  elastoplastic  solution  is  given  by 
Cruse  and  Polch  [3].  The  following  summarizes  these  developments. 

The  basic  BIE  formulation  for  a  crack  problem,  as  illustrated  in  Figure 
1 ,  is  given  as  follows 


CjiUi(P: 


+  jTj*(P,Q)ui(Q)ds  ♦ 

=  Juj*(P,Q)t.(Q)ds  + 


JT  *(P/5)u  (Q)ds 

r  -V 

*  /  0 
JU  (P,Q)t  (Q)ds 

r  J  / 
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In  (1)  the  u^,  ti  terms  are  the  boundary  displacement  and  traction  vectors  for 
the  modeled  problem.  The  kernel  functions  (or  influence  functions)  U^*, 

T^*,  are  mathematical  entities  giving  the  displacement  and  traction  that  are 
computed  on  S,  r  for  the  problem  of  an  infinite  body  loaded  at  p(x),  P(x)f  by 
a  set  of  unit  point  loads  in  each  coordinate  direction.  The  star  on  the 
kernel  functions  denotes  the  addition  to  the  point  load  solution  of  the  terms 
necessary  to  provide  for  a  traction  free  crack  at  a  specified  location  and 
orientation  in  the  geometry. 

The  use  of  a  Green's  function  for  special  geometries  is  well  developed  in 
potential  theory,  as  discussed  by  Greenberg  [4j.  In  the  current  application 
we  seek  to  obtain  fracture  mechanics  solutions  for  the  case  of  traction  free 


'Lower  case  p(x)is  an  interior  point;  upper  case  P(x) 


is  a  boundary  point. 


cracks  in  finite  planar  bodies.  The  term  with  t^Q)  for  QeT  in  (1)  is 

therefore  zero,  as  shown.  The  use  of  the  cracked  body  Green's  function 

results  in  the  traction  kernel  also  being  zero  on  the  crack,  viz. 

Tij,l(P,Q)  =  0,Q  eT,  also  as  shown  in  (1). 

Thus,  (1)  constitutes  the  constraint  equation  that  must  be  satisfied  by 

U;,  ^  on  the  uncracked  portion  of  the  surface.  This  equation  can  be  reduced 

to  solvable,  algebraic  form  through  the  use  of  suitable  approximations  to  the 

boundary  data  u^,  t^.  In  the  current  application  we  use  the  approximation  of 

piecewise  linear  interpolations  of  u^,  t^  as  developed  by  Cruse  [2]. 

The  form  of  ( 1 )  for  the  interior  displacement  provides  a  means  of  direct 

computation  of  interior  strains,  stresses,  and  stress  intensity  factors. 

Simply  stated,  the  interior  quantities  depend  on  the  totality  of  boundary  data 

for  ui ,  tj  through  integration  of  these  quantities  together  with  appropriate 

kernel  functions  for  the  cracked  plane. 

Introduction  of  inelastic  strains  (e.g.,  residual  strains  due  to  welding, 

thermal  gradient  strains,  elastoplastic  strains)  e  to  the  BIE  formulation  ft 

1  A 

results  in  a  modification  to  (1) 


A 

U 


9 


Equation  (2)  no  longer  provides  a  direct  means  for  computing  the  boundary 


.V 

A 
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V, 


data,  except  when  e  (q)  is  specified.  Thus,  for  elastoplastic  response,  an 
additional  relationship  is  needed  to  compute  the  plastic  strains  for  (2).  The 
appropriate  equation  is  the  interior  strain  distribution,  as  written  by  Cruse 
and  Polch  [3]. 

E±J(p)  =  J  S*.j(p,Q)tk(Q)ds  +  /D*iJ(p,Q)ulc(Q)ds 


+  <{>(IUm,j  +  rjim,i)£Jm(<I)dA  +  Eitmj^m(P) 


.A 


(3) 


a 


Equation  (3)  computes  the  interior  total  strain  increment  in  oerms  of  the 

boundary  data  and  the  interior  anelastic  strain  increment.  For  elastoplastic 

•  •  •  A 

solutions,  the  unknown  data  u.,  t.,  £..  are  solved  for  incrementally  and 

j  i  ij 

equations  (2),  (3)  are  coupled  on  an  iterative  basis.  The  interior  inelastic 
strains  are  modeled  as  piecewise  constant  over  AA^  area  segments  in  the 
current  study.  The  full  solution  algorithm  for  the  elastoplastic  case  is 
given  in  Figure  2.  The  yield  criterion  has  to  be  satisfied,  giving  the  amount 
of  total  strain  that  is  plastic  at  each  load  level.  The  use  of  iteration  as 
opposed  to  a  tangent  modulus  formulation  allows  us  to  precompute  all  of  the 
elastic  kernel  functions,  to  invert  one  of  these,  and  to  perform  all  of  the 
ensuing  numerics  as  matrix  multiplications. 

3-2  Stress  Intensity  Factor  Computations 

The  structure  of  (3)  has  been  investigated  by  Cruse  and  Polch  [3)  for 
interior  points  approaching  the  crack  tip.  It  was  found  that,  for  the 
elastoplastic  case  where  the  crack  tip  strains  can  exhibit  a  singularity  up  to 
l/o  (where  p  is  the  distance  from  the  crack  tip),  eq.  (3)  still  results  in 
convergent  integrals  in  (2),  (3).  However,  the  actual  strength  of  the  plastic 
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Figure  2.  Elastoplastic  Solution  Algorithm 
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strain  singularity  is  a  function  of  the  work  hardening  (see  Hutchinson  [6]) 
and  can  only  be  inferred  from  the  resulting  strain  distributions  after 
satisfying  the  flow  rule  implicit  in  (3). 

The  elastic  stress  intensity  factor  computation  for  the  BIE  formulation 
using  the  cracked  Green's  function  results  directly  from  the  elastic  version 
of  (3).  As  shown  by  Snyder  and  Cruse  [1],  the  kernels  in  the  two  boundary 
integrals  in  (3)  are  explicitly  dependent  on  the  inverse-square  root  of  the 
distance  of  p(x)  from  the  crack  tips  (+/-a).  Further,  for  the  nonsingular 
distribution  of  inelastic  strains  in  (3),  the  volumetric  kernel  has  the  same 
explicit  dependence.  Thus,  for  nonsingular,  inelastic  strains  we  obtain  the 
following  direct,  path  independent  evaluation  of  the  elastic  stress  intensity 
factors 


(KItKn)  =  -/  RiI’II(Q)u.(Q)ds  +  /  L.I’II(Q)t.(Q)d s 
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The  first  two  terms  in  (4)  are  those  previously  used  by  Snyder  and  Cruse 
[1]  and  by  Stern,  et  al.  [7].  These  are  path  independent  integrals  which 
provide  a  simple  quadrature  for  computing  Kj,  Kjj  from  any  solution  for  u^,  t^ 
on  a  path  around  the  crack,  but  excluding  the  crack. 

Equation  (4)  states  that  nonsingular,  inelastic  strains  modify  the 
elastic  Kj,  Kjj  values  in  an  equally  simple  sense  of  quadrature  when  these 
quantities  are  specified  in  the  volume  (area).  Some  examples  of  this 
quadrature  for  a  notch  plasticity  problem  will  be  discussed  below. 

In  developing  eq.  (4),  it  was  assumed  that  the  inelastic  strains  were 
nonsingular,  thus  neglecting  the  crack  tip  elastoplastic  effects.  The 
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additional  terms,  reflecting  the  higher  order  singular  behavior,  are 
represented  by  the  incremental  elastoplastic  strain  portion  of  eq.  (3) 


£ijP(p)  =  <{>  ^ i im , j *  .  e *m  ^dA 


+  EUmj  e*m  (P) 


(5) 


As  discussed  by  Cruse  and  Polch  [3],  eq.  (5)  is  dimensionally  homogeneous  for 
any  physical  singularity  in  plastic  strain  increment  as  p(x)+/-a,  but  the 
order  of  the  singularity  is  not  directly  solvable  from  (5). 

A  considerable  number  of  technologically  important  problems  to  the 
aerospace  industry  are  associated  with  the  use  of  linear  elastic  fracture 
mechanics  parameters  (i.e.,  Kj,  Kjj)  for  problems  of  limited  or  localized 
plasticity.  These  include  predicting  Kj  for  cracks  which  are  undergoing 
cyclic  plasticity  resulting  in  crack  closure  effects  on  spectrum  crack  growth, 
and  cracks  growing  in  the  plastic  zone  of  a  bolthole  subject  to  high  loading 
or  prestressing. 

The  present  research  seeks  to  shed  light  on  some  of  these  problems  by 
presenting  a  stress  intensity  factor  computation  algorithm  that  can  directly 
and  unambiguously  model  these  kinds  of  limited  plasticity  effects.  For  such 
problems,  the  solution  given  in  (4)  is  to  be  used.  The  two  boundary  data 
integrals  in  (4)  reflect  the  plastic  strain  distribution  of  the  crack  tip,  as 
well  as  other  anelastic  strains  through  the  volume  integral  in  (2).  Secondly, 
the  prior  plasticity  of  a  notch  will  affect  Kj,  Ktj  through  the  volume 
integral  of  those  nonsingular  strains  in  (4).  The  resulting  values  of  Kj,  Kjj 


are  the  plasticity  corrected  elastic  stress  intensity  factors  which  define  the 
strength  of  the  elastic  singularity  which  dominates  the  plastic  singularity. 


The  use  of  this  approach  is  obviously  limited  to  crack  tip  plastic  zones  which 
are  contained  within  the  field  of  the  elastic  singularity. 


Application  of  the  appropriate  interpolations  to  the  data  in  eqs.  (2)  and 
(3)  reduces  the  integrals  to  algebraic  form.  In  general,  the  boundary 
solution  involves  an  equal  number  of  known  (applied)  boundary  data  and  unknown 
data.  Letting  the  unknown  data  be  given  by  {*},  the  product  of  the  known  data 
and  its  coefficient  matrix  terms  by  {$},  and  the  coefficient  of  the  piecewise 
constant  plastic  strains  by  [E],  we  obtain  from  (2) 

[A]  {x}  =  {y}  h.  [E]{eP]  (6) 

Similarly,  taking  [S]  and  [D]  to  be  the  elastic  coefficient  arrays  of 
the  boundary  data,  and  [G]  to  be  the  elastic  coefficient  array  for  the 
plastic  strain,  then  eq.  (3)  becomes 

{eT}  =  t S ] { t }  ♦  t D ] {u}  +  [G] {eP}  (7) 

The  strain  superscripts  in  (6)  and  (7)  refer  to  total  (elastic  plus  plastic) 
and  plastic  values,  while  the  dots  imply  that  all  of  the  variables  are  to  be 
interpreted  in  terms  of  their  incremental  evaluation. 

The  present  BEM  algorithm  makes  use  of  the  Huber-Mises-Hencky  yield 
condition  and  associated  flow  rule.  Elastic,  perfectly  plastic  material 
response  was  modeled  throughout  this  study,  but  the  code  allows  for  a  multi- 
piecewise-linear  definition  of  a  general  stress-strain  curve. 

Following  the  approach  adopted  in  the  ADINA  code,  each  increment  in  total 
strain  is  divided  into  subincrements  (Bathe  (8]).  The  number  of  subincrements 


is  selected  to  minimise  the  error  in  the  deviatoric  stress  change  within  the 
strain  increment.  The  stress  increment  for  a  given  iterate  of  increment  in 
plastic  strain  is  then  obtained  by  an  Euler  forward  integration  of  the  flow 
rule,  according  to 
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In  eq.  (8),  the  elastoplastic  matrix  relating  the  subincrements  in  stress 
and  total  strain  is  given  for  plane  strain  by 


do  =  2G[ 6 .  6  + 


6  6.  . 
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The  current  state  of  deviatoric  stress,  S^j,  and  second  invariant  of 
deviatoric  stress,  Jg,  is  updated  within  the  subincremental  integration  of 
(8).  The  tangent  modulus,  H,  is  taken  as  the  slope  of  the  effective  stress- 
effective  plastic  strain  curve,  at  the  current  level  of  effective  stress. 

Figure  2  summari2es  the  current  iteration  algorithm  for  the  solution  of 
eqs.  (6,7).  The  coefficient  arrays  [A],  [S],  and  [D]  depend  solely  on  the 
elastic  constants  of  the  material  and  the  boundary  shape.  Thus,  they  are 
computed  once  and  stored.  The  [A]  matrix  is  inverted  prior  to  storage.  The 
interior  arrays  [E]  and  [G]  are  also  dependent  solely  on  the  elastic  constants 
and  the  interior  element  modeling.  These  are  also  computed  once  and  stored. 
Note  again  that  only  that  portion  of  the  interior  expected  to  be  inelastic 


need  be  modeled.  The  expense  of  generating  [E]  and  [G]  for  crack  problems 
dictates  that  such  limited  volumetric  modeling  be  employed. 

In  the  first  iteration  at  a  given  load  step,  the  plastic  strain  increment 
in  Figure  2  is  taken  from  the  last  load  step.  The  boundary  solution  then 
responds,  in  an  elastic  manner,  to  the  increase  in  loading.  Estimated 
interior  total  strains  are  then  calculated.  Based  on  the  new  total  strain 
increment,  the  interior  stresses  and  plastic  strains  are  computed  based  on 
satisfying  the  yield  condition  through  eq.  (8).  The  plastic  strain  increment 
is  then  updated  in  both  eqs.  (6,7)  for  a  recalculation  of  the  boundary  and 
interior  solutions. 

Absolute  convergence  of  the  strain  solution  within  each  element  is 
required  for  the  iteration  process  used.  That  is,  the  maximum  difference 
between  successive  iterates  of  the  plastic  strain  correction  term  (second  term 
on  the  right  hand  side  of  eq.  (6))  is  not  allowed  to  exceed  a  user-specified 
tolerance.  This  tolerance  has  been  selected  on  the  basis  of  its  ability  to 
relate  directly  to  the  amount  of  the  displacement  increment.  A  number  of 
numerical  experiments  with  tolerances  ranging  over  10“^  to  10"^  were  conducted 
to  test  the  sensitivity  of  the  results  to  this  value.  It  was  found  that  the 
errors  in  the  displacements,  for  a  simple  uniform  stress  test  case,  were  of 
the  order  of  the  tolerances  specified.  A  decrease  in  the  tolerance  by  an 
order  of  magnitude  generally  resulted  in  a  doubling  of  the  number  of 
iterations  required  to  achieve  convergence.  A  value  of  10“^  was  used  for  the 
notch  problem  and  a  value  of  10'^  was  used  for  the  fracture  mechanics  problem, 
in  order  to  account  for  the  higher  strain  gradient. 

3.4  Numerical  Results 

The  computer  program  has  been  verified  on  two  example  problems.  The 
first  is  a  plate  with  perforations,  previously  solved  by  Haward  and  Owen  [9] 


FT" 


■-U  \  *  T»  'A  '■  *•  T  '.V  '.’-  'A  'A'  'AH  A  'A'  AWA 1 VA  IVA'WA-A  !.WA  ■.W.V.WV.WV: 


v.V.T.V 


■f. 

V, 


,v 


V 

V 

V 


> 

A 


.y 

.y 


'A 


■r. 


using  finite  elements  and  resolved  by  Telles  [10]  using  the  BEM.  This  example 
served  the  basic  purpose  of  validating  the  current  code  and  provided  a  basis 
for  some  numerical  experimentation.  The  second  problem  is  a  fracture 
mechanics  problem  of  a  center  cracked  plate  loaded  in  tension.  The  plastic 
strain  results  are  compared  to  ADINA  results  using  a  singular  finite  element 
model . 

The  geometry  for  the  first  problem  is  shown  in  Figure  3.  Plane  strain 
conditions  are  applied  for  all  three  of  the  analyses  and  the  material  is  taken 
to  be  elastic-perfectly  plastic.  The  appropriate  constants  are  E  =  42.  x 
lO^MN/n^,  0^  =  105.  MN/m^,  v  =  0.33.  The  one  loading  condition  considered 
was  uniaxial  tension,  applied  by  prescribing  displacements  at  the  edges  of  the 
plate  section.  The  piecewise  linear  plastic  strain  BEM  mesh  of  Telles  is 
shown  in  Figure  3;  the  FEM  quadratic  isoparametric  element  mesh  used  by  Haward 
and  Owen  is  shown  in  Figure  4.  The  current  BEM  mesh,  using  piecewise  constant 
plastic  strains,  is  shown  in  Figure  5. 

Figure  6  plots  the  numerical  results  in  terms  of  the  amount  of  force 
required  versus  the  applied  displacements.  Table  1  summarizes  the  numerical 
force-displacement  data.  All  three  model  results  show  excellent  agreement, 
given  the  disparity  in  modeling  strategies.  The  predicted  limit  load  for  the 
current  study  differs  from  the  other  two  by  less  than  2%.  The  difference  is 
attributed  to  the  use  of  constant  strain  elements.  Limit  load  is  obtained 
when  the  centroidal  value  of  stress  in  the  last  ligament  element  yields,  a 
condition  that  will  occur  below  the  load  for  yielding  the  last  physical 
ligament  ahead  of  the  notch.  The  plastic  zones  are  shown  in  Figure  7  for 
various  load  levels. 

The  current  BEM  code  was  tested  for  a  range  of  load  increments  in  a 
deliberate  attempt  to  create  numerical  instability.  The  numerical  results 
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Figure  3.  Perforated  Polystyrene  Plate  and  BIE  Discretization 
Used  by  Telles  (1983) 


Results  of  the  Polystyrene  Plate 


Table  1. 


The  ultimate  strength  =  50.471  MN/m 
a,  d  defined  in  Figure  2 
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plotted  in  Figure  6  were  generated  using  an  incrementation  scheme  resulting  in 
a  single  element  yielding  at  a  time.  The  BEM  results  required  about  twenty 
iterations  per  load  step  to  fully  converge.  The  worst  case  was  one  load  step 
to  the  maximum  displacement.  The  solution  converged  in  45  iterations  and 
agreed  with  the  other  limit  load  results  within  0.2%.  The  maximum  deviation 
in  calculated  plastic  strains  was  10?  in  the  last  element  to  yield. 

The  second  example  is  a  center-cracked  plate  loaded  in  tension.  The 
total  width  of  the  plate  is  8  units,  with  a  crack  size  of  2  units.  One 
quarter  of  the  geometry  was  modeled  using  ADINA,  as  shown  in  Figures  8  and 
9.  Extremely  fine  resolution  of  the  crack  tip  elements  was  taken  in  order  to 
minimize  the  error  in  the  finite  element  solution.  The  elastic  stress 
intensity  factor  for  this  finite  element  model,  using  the  crack  opening 
displacement  at  the  quarter-point  node,  was  in  error  relative  to  handbook 
values  by  about  2%. 

The  BEM  mesh  corresponding  to  the  local  finite  element  scale  is  shown  in 
Figure  10.  The  elastic  BEM  stress  intensity  factor  results  were 
indistinguishable  from  the  handbook  results.  The  FEM/BEM  meshes  were  selected 
so  as  to  provide  about  three  decades  of  plotting  data  in  terms  of  crack  tip 
distance.  The  maximum  size  of  the  plastic  zone  was  limited  solely  for 
convenience  in  the  current  study. 

Plane  strain  conditions  were  used  for  both  of  the  models.  The  elastoplastic 
material  constants  used  were  E=2.037- lO^MN/m^,  0^3.452- 10^MN/m^,  and  v=0.27. 

The  ADINA  crack  tip  model  used  quadratic,  isoparametric  finite  elements 
with  nine  interior  strain  integration  points.  The  BEM  model  used  constant 
strain  triangles  throughout.  As  noted  above,  both  meshes  were  identical  in 
the  crack  tip  region.  The  ADINA  model  used  one  layer  of  collapsed  quadratic 
elements  adjacent  to  the  crack  tip.  This  approach  induces  a  (1/r)  type  of 
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singularity  m  the  displacement  gradient  within  this  first  xayer  of 
elements.  No  singularity  modeling  is  used  in  the  BEM  plastic  strain 
distribution . 

Loading  history  was  identical  for  both  models  and  spans  the  range  of  load 
factors  of  0.0310  to  0.2075.  A  value  of  1.0  corresponds  to  yielding  of  the 
whole  plate.  A  total  of  68  load  steps  was  used  for  both  models.  The  load 
steps  satisfy  the  conditions  of  Larsson  and  Carlsson  [11].  Simply  stated, 
these  conditions  require  that  at  most  one  element  becomes  plastic  at  each  load 
increment,  and  that  the  load  increment  should  be  smaller  than  1%  of  the  load 
corresponding  to  ^jmax=0y ‘^a-  range  of  load  factors  has  been  chosen  as 

the  range  to  go  from  yielding  the  innermost  element  to  yielding  the  outermost 
element. 

ADINA  failed  to  converge  for  the  first  step  until  the  stiffness 
reformulation  (BFGS)  procedure  was  used.  After  the  first  load  step  (requiring 
20  iterations)  the  ADINA  algorithm  with  reformulation  generally  converged  with 
five  iterations.  The  BEM  algorithm,  using  elastic  "stiffnesses,"  converged  in 
ten  to  fifty  increments  at  each  load  step  with  the  higher  numbers  occurring  at 
the  higher  load  levels.  The  total  computer  time  for  the  two  models  was 
essentially  the  same,  although  the  BEM  calculations  are  cheaper  per  load 
step.  A  higher  final  load  level  or  cyclic  loading  would  yield  a  benefit  to 
the  BEM  model,  even  though  the  current  BEM  code  is  not  yet  optimized  for  these 
calculations. 

The  crack  tip  plastic  strain  distribution  results  are  shown  in  Figure  11 
for  two  of  the  computed  load  levels.  The  data  are  taken  from  points 
distributed  near,  but  not  on,  a  line  at  an  angle  of  about  85°  to  the  plane  of 
the  crack.  This  angle  corresponds  :o  the  line  of  maximum  equivalent  elastic 
strain.  The  jaggedness  of  the  curves  is  mostly  due  to  the  use  of  triangular 


elements,  as  well  as  to  the  points  having  different  angular  locations.  The 
data  is  plotted  in  terms  of  the  centroidal  value  of  plastic  strain.  The 
innermost  row  of  finite  elements  has  three  sampling  points  radially,  which 
accounts  for  the  smaller  radius  plotted  for  these  results.  The  numerical 
results  from  ADINA  show  a  tendency  for  a  spurious  peak  in  plastic  strain 
increments  in  the  second  row  of  elements.  This  peak  is  no  doubt  induced  by 
the  lack  of  a  singularity-transition  element  in  the  current  study. 

It  is  significant  to  find  that  the  numerical  results  are  in  such  good 
agreement.  This  confirms  the  accuracy  of  the  BEM  algorithm  for  piecewise 
constant  plastic  strains.  The  BEM  results  do  not  show  the  strength  of  the 
plastic  singularity  as  strongly  in  the  first  row  of  elements  as  do  the  finite 
element  results,  with  the  imbedded  1/r  singularity  in  displacement  gradient. 
However,  both  sets  of  results  strongly  indicate  that  the  plastic  strain  for 
localized  plasticity  possesses  the  same  1/r  singularity  that  is  associated 
with  fully  developed  plasticity  for  the  case  of  zero  strain  hardening. 
Clearly,  the  presence  of  the  underlying  elastic  singularity  field  plays  an 


important  role  in  enhancing  the  modeling  accuracy  for  crack  tip  plasticity. 

Figures  12  and  13  show  the  progressive  development  of  the  plastic  zone  up 
to  the  maximum  modeled  load.  It  is  to  be  emphasized  that  the  current  study 
was  intended  to  confirm  the  accuracy  of  the  new  BEM  algorithm  for 
elastoplastic  fracture  mechanics  analysis,  and  not  study  extensive  plasticity 
response  at  the  crack  tip.  For  this  reason  the  current  results  were  not 


carried  beyond  the  load  level  shown.  There  is  no  inherent  limit  to  the  load 
level  that  can  be  modeled  with  this  BEM  algorithm. 

3-5  Crack  Extension 


The  numerical  implementation  of  the  elastoplastic  fracture  mechanics 
algorithm  in  this  section  report  focuses  on  problems  for  which  the  plasticity 
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Figure  13.  Growth  of  the  Plastic  2one  (BIE) 
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modified,  stress  intensity  factor  in  eq.  (4)  is  useful.  In  matrix  form,  this 
becomes 

Kj  n  =  <R>  {u}  +  <L>  {t}  +  <M> {eP}  (10) 

These  problems  consist  of  elastic  cracks  in  the  presence  of  plastic  strains 
due  to  welding  and  to  yielding  of  notches. 

The  first  problem  was  selected  to  validate  the  stress  intensity  factor 
algorithm  for  prior  plasticity  for  any  residual  or  thermal  strain  field.  The 
geometry  selected  is  a  simple  tension  specimen  with  the  boundary  and  internal 
mesh  shown  in  Figure  14.  The  mesh  arrangement  was  selected  solely  for  conve¬ 
nience,  as  it  is  used  as  a  portion  of  a  later  mesh.  The  specimen  was  loaded 
to  110/5  of  the  yield  stress  for  a  bilinear  material  response.  This  induced  a 
uniform  plastic  strain  throughout  the  specimen. 

The  next  step  in  the  validation  of  eq.  (10)  was  to  introduce  a  crack, 
done  along  the  bottom  of  the  mesh  as  shown.  The  residual  boundary  solution 
corresponding  to  the  residual  internal  strains  is  computed  for  the  cracked 
case  by  eq.  (6).  Next,  the  internal  strains  for  the  residual  boundary  and 
internal  variables  are  computed  from  eq.  (7).  In  the  case  of  the  test 
problem,  eq.  (6)  produced  the  uniform  displacements  compatible  with  uniform 
residual  strains;  eq.  (7)  computed  internal  strains  equal  to  the  residual 
strains . 

The  elastic  stress  intensity  factor  for  the  problem  was  then  computed 
for  the  residual  boundary  terms  computed  from  eq.  (6).  If  there  were  further 
changes  in  the  residual  strains  due  to  unloading  plasticity,  these  would 
modify  the  elastic  strain  intensity  factory  through  the  appropriate  term  in 
eq.  (10).  As  required  for  this  simple  case,  the  residual  stress  intensity 


factor  was  zero.  The  residual  values  in  the  first  and  third  terms  in  eq.  (10) 
cancel  each  other  to  within  computer  accuracy. 

Figure  15  shows  the  residual  stress  in  a  large  plate  containing  a  simu¬ 
lated  weld  bead.  The  weld  is  simulated  by  a  narrow  (b/W  =  1/400)  strip  of 
material  with  plastic  strains  equivalent  to  100  ksi  elastic  strains  in  an 
infinite  sheet.  Both  plane  strain  and  plane  stress  solutions  were  uses.  The 
plastic  strains  simulate  the  behavior  of  material  which  is  heated  to  yielding 
in  a  confined  region,  and  then  allowed  to  cool.  The  slight  accommodation  of 
the  plate  to  equilibrate  the  residual,  welding  stresses  is  seen  in  Figure  15. 

The  BIE/CRX  analysis  was  used  to  obtain  K-solutions  for  various  central 
crack  sizes  for  cracks  transverse  to  the  weld  bead.  The  only  reason  for  this 
configuration  is  to  be  able  to  compare  the  numerical  results  to  an  exact  solu¬ 
tion  for  the  residual  stresses  in  Figure  15. 

The  approach  taken  is  to  solve  eq.  (£)  to  establish  the  boundary  solution 
corresponding  to  the  residual  strains.  Two  triangles  were  used  to  integrate 
the  plastic  strains  in  eq .  (6,10);  a  single  quadr ilaterial  element  would  also 
suffice . 

The  results  of  K ( a ;  are  compared  to  the  analytical  results  using  an 
influence  function  approach  [12]  in  Figure  16.  The  agreement  is  essentially 
exact,  as  expected.  The  volume  integral  in  eq .  (10)  contributes  the  bulk  of 
the  K-solution,  as  the  boundary  displacements  associated  with  the  weld  plastic 
strains  are  very  localized. 

The  algorithm  in  eq .  (10)  is  therefore  seen  to  be  a  very  powerful  solu¬ 
tion  for  residual  plastic  strains  for  geometries  without  known  Green's  func¬ 
tions  or  influence  functions.  Further,  by  slight  reformulations  of  the  volume 
terms,  other  volumetric  strains  such  as  thermal  strains  or  body  force  strains 
can  be  analyzed  in  the  same  manner. 


Figure  15.  Simulated  Residual  Stresses  Due  to  Welding 


An  important  engineering  application  of  this  approach  is  the  calculation 
of  K(a)  for  cracks  at  notches  (holes,  fillets,  etc.)  which  have  plastic 
strains  due  to  overloads,  or  cold  working.  To  simulate  this  problem,  a  bolt 
hole  specimen  was  modeled  as  shown  in  Figure  17.  The  plate  was  loaded  in 
plane  strain  to  80?  of  the  net  section  plane  strain  yield  stress  (F^„  =  55.8 
ksi).  Figure  18  shows  the  progressive  development  of  the  plastic  strain, 
computed  by  the  BIE/CRX  program  with  zero  crack  length. 

An  elastic  crack  was  then  simulated  for  a/R  =  0.05,  0.50.  The 
solution,  as  before,  requires  an  equilibrium  ad-justment  to  account  for  the 
crack  in  the  boundary  solution  (6),  and  a  K-calculation  from  eq.  (10).  In 
this  case,  the  boundary  terms  in  (10)  contributed  a  significant  portion  of  the 
solution. 

Table  2  presents  the  elastic  and  plasticity-modified  stress  intensity 
factors  for  the  two  crack  lengths.  The  approach  used  is  to  calculate  K(a) 
from  eq.  (10)  at  no  load.  Table  2  gives  the  values  of  K(a)  at  no  load  and 
full  load  through  an  addition  of  the  elastic  result  to  the  no  load  solution. 
The  results  clearly  demonstrate  how  notch  plasticity  causes  substantial  crack 
closure  (K(a)  <  0),  thereby  reducing  the  stress  intensity  factors  at  full 
load.  Such  retardation  effects  due  to  residual  stress  have  been  previously 
modeled  with  influence  function  approaches  [13,1^]. 


/  / 


Figure  18.  Progressive  Generation  of  Plastic  Zone  at  Hole 
(Tension  Only) 


Table  2 


Stress  Intensity  Factor  Results 


o/ua 


a/R  =  0.05  a/R  =  0.5 


Elastic 

-  Max.  load  3.411  2.070 

-  Zero  load  0  0 


Plastic  (Full  Crack  Surface 
Unloading) 

-  Max.  load  0.843a  1.787a 

-  Zero  load  -2.568b  0.283b 


Note  a:  Maximum  load  values  (elastically)  =  Elastic  +  Zero  load  values 


b:  Negative  values  of  KI  imply  crack  closure  at  positive  load 


The  final  example  considers  the  effect  of  crack  tip  plasticity  on  crack 
closure.  The  selected  plane  strain  problem,  shown  in  Figure  10,  was 
previously  used  to  validate  the  BIE/CRX  algorithm  for  elasto-perfectly-plastic 
crack  tip  behavior.  Tne  plastic  strain  distribution  for  this  analysis  is 
shown  in  Figure  11,  where  a  load  level  of  one  is  net  section  yield,  in  plane 
strain . 

The  plastic  strains  in  Figure  11  clearly  indicate  the  inverse-crack-tip- 
distance  singular  behavior  expected  from  the  perfect  plasticity  solution 
[15].  As  discussed  in  [3],  such  a  plastic  strain  singularity  does  not  lead  to 
a  convergent  volume  integral  for  eq.  (4).  At  the  maximum  applied  load  (LF  i 
0.2075),  the  elastic  stress  intensity  factor  is  given  by: 


o/7ta 


with  o  =  11,579  psi  and  a  =  100  inches.  Use  of  the  full  plastic  strains  from 

the  perfect-plasticity  solution  at  this  load  level  gives  an  apparent  strength 

* 

of  the  elastic  singularity  (K^  ' 


=  1.167 


The  increase  in  Kj  using  eq.  (4)  substantially  exceeds  the  effect  of  crack 
tip  plasticity  on  crack  driving  force  that  would  be  predicted  using  a  crack 
size  equal  to  the  physical  size  plus  the  plastic  2one  size.  Denoting  the 
plastic  zone  size  by  Tp  we  obtain 
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=  0.03 


Then  the  effective  crack  tip  stress  intensity  factor  would  be  estimated  as 


=  1.015 


Thus,  we  conclude  that  the  singularity  contribution  to  eq.  (4)  is,  in  fact, 
unbounded  and  the  result  in  (12)  is  not  meaningful. 

The  unloaded  solution  to  the  same  crack  problem  involves  reversal  of  the 
crack  tip  plasticity.  Figure  19  compares  the  loaded  and  unloaded  transverse 

P 

plastic  strain  (e^)  distributions  along  the  line  of  maximum  equivalent  stress 
The  unloaded  plastic  strain  distribution  is  much  less  singular  and  we  compute 
at  zero  load 


=  0.11 


as  compared  to  the  value  0.015  estimated  in  (14).  This  large  result  still 
indicates  that  the  residual  strain  singularity  is  still  too  strong  to  be 
neglected  in  using  eq .  (4). 


Figure  19.  Loaded  and  Unloaded  Transverse  Plastic  Strain 


While  the  magnitude  of  (15)  is  not  meaningful,  eq.  (4)  can  be  used  to 
indicate  the  relative  effect  of  the  residual  strains  on  retardation  of  crack 


growth.  The  approach  is  the  same  as  used  in  the  notch  problem  discussed 


above.  That  is,  the  crack  is  extended  elastically  into  the  residual  strain 

# 

# 

field,  and  KT  is  recomputed.  The  result,  in  Figure  20,  is  given  as  r; — . 

KI 


The  plot  shows  almost  an  instantaneous  change  of  K.  at  no  load  to  a 
negative  value.  This  result  confirms  the  expectation  of  a  retardation  effect 


of  crack  tip  plasticity  on  subsequent  cyclic  stress  intensity  factor.  The 


peak  effect  is  at  —  =  0.27,  corresponding  roughly  to  some  estimates  of  the 
rp  r 

plane  strain  plastic  zone  size  the  extent  of  the  retardation  zone 


ls-=3. 


Thus,  we  conclude  that  the  volume  integral  in  eq.  (4)  is  substantially 
nonconvergent  for  the  crack  problem.  The  correction  would  be  to  delete  from 
the  plastic  strains  those  due  to  the  singularity  alone.  It  is  expected  that 
the  remaining  plastic  strains  can  be  used  to  compute  an  effective  stress 
intensity  factor  which  accounts  for  contained  plasticity  effects. 


4.0  GENERAL  SOLUTION  PROCEDURE  FOR  FRACTURE  MECHANICS 
WEIGHT  FUNCTION  EVALUATION 


4 . 1  Weight  Functions 

The  weight  function  method  is  based  on  Rice's  [16]  interpretation  of 
Bueckner's  [17]  original  paper.  The  weight  function  for  a  crack  problem  is 
generally  taken  to  be  the  normalized  rate  of  change  of  surface  displacements 
with  respect  to  crack  size  for  a  reference  state  of  loading.  As  shown  by  Rice 
[16],  this  weight  function  acts  as  a  Green's  function  for  the  crack  problem. 
That  is,  the  solution  to  any  fracture  mechanics  problem  for  the  same  geometry 
but  different  loading  conditions  can  be  obtained  from  the  weight  function  for 
the  reference  set  of  loading  conditions.  The  manner  in  which  this  is  done 
will  be  reviewed  in  greater  detail  later  in  this  paper,  but  the  process 
involves  an  integration  of  the  uncracked  stress  field  times  the  weight 
function  to  arrive  at  the  crack  tip  stress  intensity  factor  for  those  imposed 
stresses. 

The  singular  advantage  of  the  weight  function  method  is  efficiency  of 
computation  of  the  crack  tip  stress  intensity  factor  for  a  variety  of  crack 
sizes  and  loading  conditions.  Crack  size  effects  such  as  finite  width  effects 
on  stress  intensity  factor,  or  size  effects  where  crack  size  changes  the 
applied  loads  (stiffness  effects),  need  to  be  included  in  the  weight  function. 
The  new  method  reported  herein  addresses  the  computational  problem  of  gener¬ 
ating  weight  functions  in  a  direct  and  efficient  manner,  while  providing  for 
the  first  time  a  general  method  for  the  mixed-mode  problem. 

4.2  Numerical  Methods  for  Evaluating  Weight  Functions 

The  weight  function  method  discussed  herein  excludes  the  related,  Green's 
function  method  for  stress  intensity  factor  evaluation.  In  the  Green's 
function  method,  many  of  which  are  given  in  Rooke  and  Cartwright  [18],  the 


stress  intensity  factor  for  the  reference  problem  is  given  in  generalized 
terms  as  the  response  to  a  point  load  on  the  crack  surface.  The  point  load 
solution  is  then  integrated  as  a  weighted  integral  of  the  applied  tractions. 

While  the  Green's  function  and  weight  function  methodologies  are  closely 
related,  the  computational  approaches  for  the  two  methods  are  distinct.  The 
current  paper  focuses  on  the  numerical  evaluation  of  surface  displacement 
derivative  evaluation  to  establish  the  weight  function.  The  approach  is  based 
on  the  boundary  element  method  for  two-dimensional  fracture  mechanics,  as 
first  developed  by  Snyder  and  Cruse  [1],  enhanced  and  corrected  by  Cruse  [2], 
and  most  recently  used  for  elastoplastic  fracture  mechanics  modeling  by  Cruse 
and  Polch  [3].  In  these  reports  the  general  BEM  is  implemented  with  augmented 
boundary  integral  equation  kernels  which  explicitly  include  the  presence  of  a 
stress  free  crack. 

Most  numerical  methods  for  the  development  of  weight  functions  are  based 
on  the  boundary  integral  equation,  Besuner  [19],  or  finite  element  methods, 
Parks  [20],  The  boundary  integral  equation  solutions  (now  referred  to 
generally  as  boundary  element  methods)  have  been  based  upon  numerical  differ¬ 
entiation  of  the  numerical  results.  That  is,  the  crack  surface  displacement 
numerical  solutions  are  obtained  for  two  slightly  different  crack  lengths. 

The  finite  element  method  was  modified  by  Parks  [20]  to  include  in  the  virtual 
work  principle  an  explicit  derivative  with  respect  to  crack  size.  This 
permitted  the  FEM-based  approach  to  be  more  efficient  than  numerical  differ¬ 
entiation.  This  approach  has  been  exploited  to  a  great  extent  by  Sha  and  Yang 
[21]. 

A  second  approach  to  the  numerical  problem  is  that  proposed  originally  by 
Paris,  McMeeking  and  Tada  [22].  In  this  method,  the  weight  function  is 
computed  for  the  problem  of  cracked  body  subject  to  the  elastic  singularity 


tractions  on  a  small  circular  path  surrounding  the  crack  tip.  While  most  of 
the  results  were  obtained  using  the  FEM,  Cartwright  and  Rooke  [23]  used  the 
BEM  to  obtain  very  satisfying  numerical  results.  Others,  e.g.,  Grandt  [24] 
and  Petroski  and  Achenbach  [25],  have  developed  very  efficient  means  of 
estimating  weight  functions  for  a  limited  number  of  specific  geometries. 

4 . 3  Formulation 

A  general  purpose  numerical  evaluation  of  crack  surface  weight  functions 
has  been  developed  based  on  a  novel  use  of  a  specialized  boundary  element 
formulation  of  the  two-dimensional  fracture  mechanics  problem.  The  following 
sections  will  present  two-dimensional  brief  definition  of  crack  surface  weight 
functions,  their  use,  and  some  of  their  potential  limitations,  followed  by  a 
review  of  the  new  method.  The  basic  references  for  these  discussions  are  Rice 
[16]  for  the  weight  functions,  and  Cruse  [2]  for  the  boundary  element 
method.  The  general  weight  function  approach  of  Bortman  and  Banks-Sills  [26] 
will  be  followed. 

4.3.1  Crack  Surface  Weight  Function 

Consider  two  solutions  for  the  specified  geometry.  We  will  call 
solution  state  -(2)  as  the  unknown  solution,  and  solution  state  -(1)  as  the 
reference  state.  Both  solutions  consist  of  stress  and  strain  variables  on  the 
interior  of  the  body,  and  tractions  and  displacements  on  the  surface.  The 
reciprocal  work  done  by  the  stresses  of  solution  -(2)  on  the  strains  of  the 
reference  state  is  given  by* 
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•Superscripts  should  not  be  confused  with  powers  in  this  section. 


The  boundary  integral  consists  of  the  tractions  t.  and  the  displacements  uT" 

2 

for  the  two  solutions;  f  is  the  body  force  term. 

Consider  next  the  virtual  extension  of  the  crack  for  the 

reciprocal  problem.  The  tractions  that  are  released  ahead  of  the  crack  in 

(16)  are  singular  with  the  usual  inverse  square-root  behavior.  The 

displacement  term  in  (16)  is  proportional  to  the  square  root  of  crack  tip 

distance.  A  reciprocal  energy  release  rate  can  be  computed  in  terms  of  the 

mixed-mode  stress  intensity  factors  (Kj ,  )  for  the  two  solutions.  Letting 

2 

H  =  E/( 1  -  v  )  for  plane  strain  and  H  =  E  for  plane  stress,  we  obtain 
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The  boundary  has  been  divided  into  that  part  for  specified  tractions  (S  )  and 
that  for  specified  displacements  (Su).  The  body  force  term  in  (17)  has  been 
dropped  only  for  simplicity  at  this  point.  Since  state  -(2)  is  general,  we 
can  compute  equation  (17)  for  another  general  (and  independent)  solution.  Th 
reciprocal  strain  energy  release  rate  in  (17)  then  gives 
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Combining  equations  (17)  and  (18)  we  can  solve  for  the  stress  intensity 
factors  for  the  arbitrary  problems  in  terms  of  the  solution  of  the  reference 


problem 
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( 19b) 


where 


(K)2  =  K2  -  K2j  *  0 


In  the  case  of  symmetric  loading,  equation  (20)  reduces  to 
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This  is  the  same  form  as  developed  by  Rice  [16]  except  for  the  introduction  of 
mixed  boundary  conditions.  The  solution  approach  using  weight  functions  is 
normally  given  for  =  0  such  that 
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where  h^S)  is  the  weight  function.  For  the  mixed  boundary  value  problem,  no 
simple  weight  function  can  be  written  down  and  the  full  form  of  equation  (21) 
has  to  be  used. 
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Thus,  for  the  general  weight  function  method,  one  needs  an 
efficient  means  for  solving  the  reference  geometry  fracture  mechanics  problem, 
subject  to  any  fixed  displacement  or  compliant  boundary  conditions,  to  obtain 
the  rate  of  change  in  boundary  displacements  and  tractions,  and  the  crack  tip 
stress  intensity  factors.  The  fracture  mechanics  BEM  provides  the  most  direct 
and  efficient  means  for  providing  this  data  in  the  appropriate  manner. 

4.3.2  Boundary  Integral  Equation  Formulation 

The  boundary  element  formulation  of  the  two-dimensional  fracture 
mechanics  problem  is  restated  for  completeness  in  a  shortened  notation 


u/2  +  j  T*u  dS  =  f  U*t  dS 

s  s 


where  the  physical  variables  are  the  boundary  displacement  and  traction 
vectors,  u,  t.  The  kernel  functions  U,  T  in  (23)  are  the  displacement  and 
traction  (on  S)  solutions  to  the  elasticity  problem  for  the  infinite  plane 
subject  to  a  point  force  loading;  these  functions  are  known  as  fundamental 
solutions  in  the  BEM.  The  asterisk  on  each  of  the  fundamental  solutions 
denotes  further  that  they  represent  the  solution  for  the  point  load  in  an 
infinite  plane,  containing  a  single,  traction  free  crack  of  length  2a.  The 
surface  S  in  (24)  does  not  include  the  crack  surface,  as  these  boundary  condi¬ 
tions  are  automatically  satisfied  by  the  special  Green's  function  or  funda¬ 
mental  solution  used  in  this  formulation. 

The  weight  function  method  requires  the  solution  in  terms  of  the 
rate  of  change  of  the  boundary  conditions,  as  a  function  of  crack  length. 
Differentiation  of  (23)  with  respect  to  crack  length  is  completely  straight¬ 
forward,  due  to  the  explicit  dependence  on  crack  length  of  the  kernels  of  the 
integral  operators.  Thus,  we  obtain  the  following  boundary  identity 
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=  [A]{x]  =  [dB/da] {y} 


(25) 
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The  boundary  has  been  denoted  as  S  =  S, ,  +  Sf  to  denote  the  notions  of  both 
displacement  and  traction  boundary  conditions  on  portions  of  the  surface.  In 
general,  the  mixed  boundary  conditions  involve  somewhat  more  complexity  than 
is  denoted  by  this  notation,  but  the  terms  in  (24)  convey  the  essential  notion 
of  the  algorithm. 

Letting  the  vector  x  correspond  to  the  unknown  boundary  conditions 
(both  traction  and  displacement  components),  and  y  the  known  boundary  condi¬ 
tions,  a  symbolic  form  of  (24)  may  be  written 


Equation  (25)  is  obtained  from  (24)  by  the  imposition  of  a  boundary  interpola¬ 
tion  system.  In  the  current  application,  the  boundary  data  is  assumed  to  vary 
in  a  piecewise  linear  fashion,  with  x,  y  representing  nodal  variables. 

Equation  (25)  is  very  similar  to  the  discrete  form  of  the  BIE  (23) 
which  is  solved  to  obtain  the  unknown  boundary  conditions  for  the  reference 
problem.  Specifically,  A  is  identically  the  same  as  for  the  boundary  value 
problem.  The  right-hand-side  of  (25)  is  made  up  of  the  derivative  of  the  B 
terms  with  respect  to  the  crack  length,  and  y  is  the  totality  of  boundary  data 
from  the  solution  of  the  reference  problem,  equation  (23). 

Thus,  the  algorithm  for  the  solution  of  (25)  consists  of 

1.  setting  up  the  discrete  form  of  (23)  to  obtain  A, 

2.  storing  A  for  later  use, 

3.  solving  for  the  unknown  boundary  data, 
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4.  computing  the  elements  of  [dB/da], 

5.  solving  (25)  using  the  entirety  of  the  boundary  conditions 
The  solution  to  (25)  consists  of  the  change  in  all  non-specif ied 

boundary  conditions,  except  along  the  crack  surface.  This  is  so  because  the 
special  Green's  function  used  for  the  fracture  mechanics  solution  automati¬ 
cally  satisfies  the  crack  surface  conditions,  and  the  crack  surface  is  elimi¬ 
nated  from  (23).  The  complete  solution  therefore  requires  an  additional  step 
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before  the  full  weight  function  can  be  defined. 

The  BEM  is  based  on  the  so-called  Somigliana  identity  for  the 
interior  solution  variable,  which  in  this  case  is  the  displacement  variable. 
The  derivative  of  the  interior  displacements  corresponding  to  the  boundary 
solution  co  (23)  is  given  by 
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(26) 


Equation  (26)  is  equally  true  for  interior  points  or  crack  surface  points, 
although  care  must  be  exercised  in  the  evaluation  of  the  kernel  functions  for 
crack  surface  points.  The  numerical  solution  of  (26)  completes  the  weight 
function  description,  in  a  formal  way.  That  is,  by  specifying  a  suitable  num¬ 
ber  of  crack  surface  locations,  the  analyst  is  given  a  complete  description  of 
the  rate  of  boundary  condition  change  with  respect  to  crack  length,  as  needed 
in  equation  (21).  However,  the  term  from  (26)  for  the  point  at  the  crack  tip 
is  singular  and  requires  special  treatment,  as  discussed  in  the  next  section. 

4.3.3  Crack  Tip  Stress  Intensity  Factor 

One  of  the  very  important  features  of  the  fracture  mechanics  BIE 
solution  is  direct  access  to  a  very  accurate  and  generic  algorithm  for  eval¬ 
uating  the  crack  tip  stress  intensity  factors  for  mixed  mode  response.  This 
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algorithm  has  been  previously  discussed  and  requires  only  that  the  following 
integrals  be  evaluated  for  some  path  in  the  cracked  body 


% II 


=  - 1  Bi,nu  dS*|Li,nt 


dS 


(27) 


The  complete  traction  and  displacement  solution  needs  to  be  specified  on  some 
path,  and  the  original  surface  (excluding  the  crack  surface)  is  usually 
selected.  However,  as  discussed  in  Cruse  [27],  the  algorithm  is  independent 
of  path  and  has  successfully  been  used  as  a  post-processing  algorithm  even  for 
finite  element  models. 

The  weight  function  algorithm  also  involves  the  direct  computation 
of  a  crack  tip  singularity.  Taking  the  free  term  in  (26)  to  be  evaluated  near 
the  crack  tip,  it  has  been  found  that  the  kernels  contain  an  explicit,  square- 
root  singularity.  As  in  the  development  of  (27)  we  take  the  integral  identity 
(26)  for  points  near  the  crack  tips  (either  end  of  the  crack),  multiply  by  the 
square  root  of  crack  tip  distance  and  proceed,  in  the  limit,  to  the  crack  tip. 
A  non-singular  form  of  (26)  is  then  obtained 
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These  terms  are  directly  related  to  the  stress  intensity  factors  found  in 
(27),  but  are  derived  from  a  completely  independent  basis.  Thus,  they  provide 
a  validation  of  equation  (25)  by  comparison  to  the  results  from  equation  (27). 
Having  the  explicit,  limiting  singularity  strength  provides  the  analyst  with  a 
complete  description  of  the  weight  function  for  the  reference  problem. 


4.4  Applications 


The  applications  given  in  this  paper  are  limited  to  validation  examples 
for  the  computational  algorithm  that  has  been  developed.  No  effort  to  compute 
stress  intensity  factors  from  the  BEM-developed  weight  functions  will  be 
given,  as  these  are  covered  adequately  in  the  references  given.  The  following 
sections  will  present  a  brief  description  of  the  computational  implementation 
of  the  BEM-weight  function  algorithm,  two  numerical  examples  of  interest,  and 
some  conclusions  based  on  the  work  done  to  date. 

4.4.1  Programming 

The  BIE  (23)  is  first  formulated,  making  direct  use  of  the  defined 
boundary  conditions  to  assemble  the  matrix  for  the  unknown  boundary  data.  The 
right  hand  side  vector  is  assembled  from  the  product  of  the  kernel  functions 
and  the  specified  boundary  data.  The  coefficient  matrix  is  stored  on  tape  for 
later  recall,  and  the  system  of  equations  is  solved  by  reduction.  The  full 
set  of  boundary  data  is  then  completed  and  the  coefficient  matrix  recalled. 

The  new  right  hand  side  vector  is  formed  from  the  matrices  resulting  from 
differentiation  of  the  standard  kernels,  with  respect  to  crack  length,  times 
the  full  set  of  boundary  data  for  the  reference  problem.  Re-reduction  of  the 
same  matrix  then  results  in  the  set  of  boundary  data  corresponding  to  rates  of 
change  with  respect  to  crack  length.  Mixed  boundary  conditions  and  mixed-mode 
cracking  problems  are  handled  directly. 

The  next  step  is  to  compute  the  crack  tip  terms  at  one  or  both  of 
the  crack  tips  and  the  rate  of  change  of  crack  surface  displacements  with 
respect  to  crack  length.  The  same  matrices  as  used  for  the  boundary  terms 
above  are  computed  for  the  points  located  on  the  crack  surface.  These  points 
are  inDut  as  percentages  of  crack  length,  and  are  taken  to  be  on  the  top  sur¬ 
face  of  the  crack  for  symmetric  problems. 
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4.4.2  Numerical  Examples 

Two  example  reference  problem  calculations  are  presented  which 
illustrate  the  essential  numerical  capabilities  of  the  weight  function 
calculation  procedure.  The  first  is  more  in  the  way  of  a  validation  example, 
as  the  problem  considered  is  a  simple  square  plate  with  a  central  crack, 

Figure  21.  The  plate  is  loaded  in  simple  tension  transverse  to  the  crack  as 
shown;  this  loading  results  in  Mode  I  response  of  the  crack.  The  crack  sizes 
range  from  a/W  =  0.01  to  0.5. 

The  smallest  crack  size  is  essentially  equal  to  the  infinite  plate 
problem.  The  computer  code  calculated  stress  intensity  factors  for  each  case; 
for  the  short  crack  case  the  accuracy  to  the  infinite  plate  result  was  to  five 
significant  figures.  Figure  22  plots  the  normalized  crack  opening  displace¬ 
ment  derivative  results  for  the  four  cases.  It  is  seen  that  in  all  cases  the 
resulting  distribution  is  quite  smooth  over  most  of  the  crack;  the  finite 
width  effect  is  seen  for  the  cases  of  a/w  >  0.1.  The  normalized  crack  tip 
singular  behavior  is  also  seen  in  Figure  22  to  be  essentially  identical  for 
each  of  the  cases.  The  stress  intensity  factor  computed  for  the  crack  opening 
displacements,  using  equation  (28),  was  essentially  identical  to  the 
previously  exploited  algorithm  based  on  internal  stresses,  equation  (27). 

Figure  23  shows  the  BEM  mesh  for  the  second  reference  problem 
considered  in  this  study.  The  problem  is  a  plate  with  a  central  hole  and  an 
edge  crack  from  the  hole.  The  mesh  was  established  in  a  manner  that  repre¬ 
sented  mixed  boundary  conditions,  as  the  left-hand  side  is  taken  as  a  plane  of 
simulated  symmetry. 

The  crack  is  taken  to  be  located  at  the  horizontal  symmetry  line, 
extending  from  the  edge  of  the  crack.  It  is  to  be  re-emphasized  that  the  BEM 
algorithm  being  used  does  not  model  the  surface  of  the  crack,  as  this  surface 
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is  explicitly  and  exactly  accounted  for  in  the  formulation  of  the  integral 
equations.  The  plate  is  taken  to  be  sixteen  units  long,  with  a  width  of  eight 
units  and  a  hole  radius  of  two  units. 

Figure  24  plots  the  computed  stress  intensity  factors  for  five 
crack  sizes  (a/R  =  1.5,  1,  0.5,  0.25,  and  0.05).  The  stress  concentration 
factor  for  the  hole  is  computed  to  be  3-59  versus  the  value  of  3.54  from 
Peterson  [15].  The  stress  intensity  factor  for  an  infinitesimally  short  edge 
crack  at  the  hole  is  then  given  by 

Kj  =  1.12  Kt  o  /^a  (29) 

The  value  of  normalized  stress  intensity  factor  in  Figure  24  for 
a/R  =  0,  therefore,  is  3.98.  The  stress  intensity  factor  is  seen  to  decrease 
with  increasing  crack  size  for  short  crack  lengths  due  to  the  effect  of  the 
stress  gradient.  For  longer  crack  lengths,  the  values  are  seen  to  begin  to 
rise,  as  would  be  expected  due  to  finite  width  effect.  The  second  curve  in 
Figure  24  simply  normalizes  the  stress  intensity  factor  such  that  the  value  at 
a/R  =  0  is  zero. 

Figures  25  and  26  present  the  numerical  results  for  the  normalized 
weight  function  for  the  Mode  I  reference  problem.  The  results  in  Figure  25 
are  essentially  of  the  same  shape  as  in  Figure  22,  although  the  zero  intercept 
is  elevated  due  to  the  free  edge  effect.  Also  the  order  of  the  curves  is 
affected  by  a/R.  Figure  26  normalizes  these  results  by  the  numerical  results 
for  the  center  cracked  panel  (a/W  =  0).  The  case  for  a/R  =  0.05  clearly  shows 
the  influence  of  the  free  edge  effort  for  short  cracks  while  the  others  begin 
to  look  more  like  the  results  for  a  center  cracked  panel  of  crack  length  equal 


Figure  24.  Stress  Intensity  Factions:  Tension  Plate  with  Hole 


a/4  =  0.05 


5.0  BURIED  CRACK  ANALYSIS  WITH  AN  ADVANCED  TRACTION  BIE  ALGORITHM 


V  % 


5 . 1  Introduction 

The  algorithm  for  analysis  of  buried  cracks  of  arbitrary  shape,  presented 
in  this  section,  is  a  part  of  a  bigger  code.  The  ultimate  goal  of  the 
principal  code  is  a  solution  of  an  arbitrary  surface  crack  problem.  The 
alternating  boundary  traction  method  (Nishioka  and  Atluri  [28])  is  used  to 
this  end.  Both  the  literature  and  our  experience  show  the  efficiency  and  good 
convergence  properties  of  the  alternating  method.  The  accuracy  of  the 
solution  of  a  buried  crack  problem,  being  the  main  part  of  the  alternating 
method  algorithm,  plays  the  most  important  role  in  the  accuracy  of  the  overall 
results. 

The  problem  of  a  buried  crack  with  arbitrary  geometry  has  not  been  solved 
satisfactorily.  Nishioka  and  Atluri  [28]  used  a  solution  based  on  Jacobi 
potential  functions  for  the  elliptical  planar  geometry  and  polynomial 
loading.  Weaver  [29]  used  a  dislocation  model  to  solve  the  problem  for  the 
case  of  the  rectangular  crack.  Using  the  same  dislocation  method,  Bui  [30] 
presented  a  solution  for  an  arbitrary  geometry,  but  with  poor  accuracy.  The 
accuracy  of  Bui's  method  has  been  improved  by  Putot  [31],  albeit  with  certain 
artifices.  The  current  work,  reported  herein,  removes  these  limitations  and 
deficiencies  by  a  successful  combination  of  quite  well  known  numerical 
techniques . 

5.2  Mathematical  Statement  of  the  Problem 

The  problem  of  a  loaded  crack  in  an  infinite  three-dimensional  body  can 

be  mathematically  described  by  a  set  of  integral  equations  of  the  first  kind 

£ 

for  relative  crack  opening  displacements  in  terms  of  crack  surface 
tractions  a.-. 
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where 


r2(P,Q)  =  I  (x,(P)  -  x ■ ( Q ) ) 2  ;  v.  =  v.  ( 
i=1  1  1  1>a  i*0 


<r>  denotes  integration  over  r  in  the  Cauchy  Principal  Value  sense. 


This  well-known  set  of  equations  (Weaver  [29]),  Bui  [30],  Cruse  [27])is 


herein  called  the  Traction  Boundary  Integral  Equation.  The  Greek  subscripts 


apply  to  the  in-plane  directions  defined  by  the  normal  direction,  X^.  This 


set  of  equations  is  derived  from  the  Somigliana  identity  for  the  internal 


stresses  in  a  body  with  the  crack  (Cruse  [27]).  An  important  feature  of  this 


representation  of  the  physical  problem  is  that  it  decreases  the  dimensionality 


of  the  problem  from  three  to  two,  the  equations  extending  only  over  the  two- 


dimensional  crack  surface. 


Vj^  =  (upper  surface) 


u^  (lower  surface) 
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An  important  aspect  of  the  set  of  equations  (30)  is  the  presence  of 
principal  value  integrals.  The  evaluation  of  principal  value  integrals 
usually  involves  special  treatment  of  exclusion  circles  for  singularity. 

Proper  handling  of  this  practical  problem  requires  continuity  of  displacement 
gradients  (crack  surface  strains)  between  the  elements.  If  the  continuity 
condition  is  not  satisfied,  certain  artifices  have  to  be  used  to  contain  the 
error  (Putot  [31]). 

5.3  Current  Numerical  Method 

In  the  current  numerical  method,  we  use  a  finite  element  interpolation 
scheme  for  the  unknown  displacement  discontinuities  at  the  crack  surface,  and 
their  derivatives.  In  the  following,  we  will  use  the  term  "displacements"  for 
"displacement  discontinuities"  for  brevity.  The  finite  element  interpolations 
produce  undesirable  discontinuities  of  displacement  derivatives  at  boundaries 
of  the  elements  unless  higher  order  derivatives  are  used  as  nodal  quantities. 

A  special  interpolation  procedure  is  used  to  remedy  the  problem.  We 
model  the  displacement  field  on  the  crack  surface  using  quadratic  8-noded 
isoparametric  elements,  and  the  displacement  gradient  field  using  linear,  4- 
noded  isoparametric  elements.  The  relationship  between  these  two 
interpolations  is  established  by  the  following  schemes. 

The  interpolation  of  the  displacement  discontinuity  field  on  a  crack  is 
given  by  the  relation 

v.(P)  =  Au.(P)  =  Nk(P)Au.k=Nk(P)vk  (3D 

where 

k  =  1  ,n 

P  =  point  on  the  crack  surface 
n  =  number  of  nodes  of  8-noded  element  meshes 


We  will  call  the  displacement  interpolation  mesh  using  the  8-noded  element  the 


N-mesh  in  the  following  discussion.  Differentiation  of  with  respect  to  the 
global  variable  x^  will  provide  the  in-plane  displacement  gradient  field 
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(32) 


This  field  will  be,  in  general,  discontinuous  at  the  element  boundaries  (as 
the  stresses  and  strains  are  in  standard  finite  elements). 

To  find  a  continuous  displacement  gradient  field  corresponding  to  the 
original  displacement  field,  we  will  create  a  new  mesh  of  linear  4-noded 
isoparametric  elements,  referred  to  as  the  M-mesh.  This  new  mesh  will 
interpolate  displacement  gradients  on  crack  surface  according  to  the  following 


=  MJ(P)v.  J 
l.a 


(33) 


where 

j  =  1  ,m 

m  =  number  of  nodes  of  4-noded  element  mesh 
The  nodes  of  the  M-mesh  will  coincide  with  those  of  the  N-mesh  and  with  the 
centers  of  8-noded  elements  in  a  manner  shown  in  Figure  27. 

We  will  find  the  nodal  values  of  displacement  gradients  by  minimizing  the 
difference  between  the  two  interpolations,  (32)  and  (33).  The  difference 
between  the  interpolations  on  the  whole  crack  surface  (indices  k,j  running 
through  all  nodes  of  the  crack  model)  can  be  written  in  shorter  form 
as  E.  where 

la 
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k  r  1,n 
j  =  1  ,m 
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Figure  27.  Interpolation  Meshes  on  Crack  Surface 


V  ■-  A.  v. 


r^.TT7.^^.  r.  r,  ■ 


69 


fr 


5 


A 
•<■  ‘ 


Using  the  Galerkin  weighted  residual  process,  with  the  same  weighting 
functions  as  the  displacement  gradient  interpolation  shape  functions,  we 
obtain  the  following 


f  E-MrdS  =  0 
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r  =  1  ,m 


(35; 
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Substitution  of  equation  (5)  results  in 
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The  matrix  representation  of  equation  (35b)  takes  on  the  following  form 


[MM|{v.ia}  =  [NM]  fv.} 


where 


(36) 
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[mm]  =  square  matrix  of  dimensions  6m  x  6m  with  MM.  =  J  M*'MrdS 

Jr  S 

3Nk  r 

[ MM ]  =  rectangular  matrix  of  dimension  6m  x  3n  with  NM  =  /  — r—  MrdS 
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(v.  }  =  column  vector  of  nodal  displacement  gradients 

’  (dimension  6m  x  1 ) 


{v.}  =  column  vector  of  nodal  displacements 
(dimension  3n  x  1 ) 
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Matrix  MM  has  the  same  structure  as  the  finite  element  consistent  mass  matrix; 
clearly,  MM  is  positive  definite.  Nodal  values  of  the  displacement  gradients 
are  then  calculated  by 


{v^}  =  [MM-1]  [NM]  {v.} 


This  process  very  closely  resembles  global  stress  smoothing  used  sometimes  in 
the  finite  element  method,  Hinton  and  Campbell  [32].  Since  MM  is  positive 
definite  and  has  a  strongly  dominant  diagonal,  its  inversion  does  not  present 
any  problems.  In  short,  the  global  relationship  between  nodal  displacement 
gradients  and  displacements  can  be  written  as 


{vi,a}  =  1^1  (vi^ 

(38) 

where  [ MNM]6m  x  3n  --  [MM-1]  [NM] 

The  numerical  solution  of  the  traction  BIE,  equation  (30),  uses  the 
collocation  method  to  form  the  equivalent  system  of  algebraic  equations  for 
the  unknowns.  Since  the  displacement  mesh  contains  n  nodes  with  three  unknown 
components  of  displacements  each,  we  take  these  nodal  locations  as  collocation 
points. 

The  continuous  representation  of  displacement  gradients  given  by  (38) 
allows  for  very  easy  and  natural  treatment  of  the  principal  value  integrals 
appearing  in  (30).  It  can  be  easily  proven  that  the  basic  component  of  all 
integrals  in  (30)  vanishes  as  the  radius  of  an  exclusion  circle  around  the 
source  point  P  goes  down  to  zero,  viz. 


I;  (P)  =  Lim  J  V.  (Q) 
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Discretization  of  the  set  of  equations  (30)  with  the  interpolation  scheme 
(33)  at  n  collocation  points  gives  the  following  system  of  equations 


[»»]  •  =  lo13} 


i  =  1,3 

a  =  1,2 

1  =  1,3 


where 


VV  -  3n  x  6m  array  of  coefficients 


{Vi  ^}  -  6m  x  1  column  vector  of  unknown  displacement  gradients 


[013}  -  3n  x  1  column  vector  of  known  surface  stresses  (tractions) 


This  underdetermined  system  cannot  be  solved  directly  for  the  displacement 
gradients.  To  reduce  the  number  of  unknowns,  we  use  relation  (38),  thus 
changing  the  unknowns  to  displacements.  We  obtain 


[w]  .  [MNM]  •  {v.}  =  {a13} 


[VK]  •  {v.}  =  {a  } 
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is  the  final  array  of  coefficients.  System  (42)  has  3n  displacement 


components  as  unknowns  and  3n  equations  formulated  at  n  nodes  of  the 


displacement  mesh.  It  can  be  solved  by  any  convenient  numerical  technique 


after  application  of  displacement  boundary  conditions  (43)  at  the  contour  of 


crack  surface 


v  =  0  (i  =  1,3;k  =  node  numbers) 


The  structure  of  the  system  of  equation  (42)  shows  that  (for  the  current 


formulation  for  plane  cracks)  the  problem  may  be  separated  into  two  separate 


systems  of  equations.  Crack  opening  displacements  are  not  coupled  through  the 


equations  with  in-plane  displacements ;  thus  the  system  effectively  splits  into 


two  separate  systems  of  equations 


[VK  ]  •  {v.}  =  {o.J 

Ayj  i  j_j  i  13J 


i  =  1,2 

1  =  1,2 


(43a) 


Kl  •  {v3i  =  t°33i 


(43b) 


This  property  allows  for  greater  efficiency  of  numerical  solution-systems  of 


equations  for  separate  sub-problems  are  smaller  than  the  full  system.  The 


decoupling  of  the  system  is  used  in  the  computer  implementation  of  the 


method.  After  solution  of  system  (42)  the  displacement  gradients  at  m  nodes 


of  the  linear  element  mesh  can  be  obtained  using  equation  (38). 


Stresses  at  any  location  of  the  infinite  body  due  to  crack  loading 


(except  for  the  crack  surface  itself)  can  then  be  obtained  using  formula  (44) 


°ia(p)  =  -  J  slia(p,Q)v1(Q)dS(Q) 
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The  kernel  S 


is  given  in  Cruse  (333.  To  calculate  integral  (44),  we  use 


lia 

the  displacement  interpolation  (31)  on  the  M-mesh.  Integral  (44)  presents  no 
particular  problems,  since  it  is  a  proper  one.  Higher  order  numerical 
integration  is  only  necessary  when  the  source  point  p  is  located  close  to  the 
crack  surface. 

Stress  intensity  factors  at  nodes  of  crack  front  are  calculated  using  the 
formulas  of  Bui  [30]: 


VIE 


1  8<1-v2) 


(45) 


LII 


VHE 

8(1-V 


(46) 


VIIIE 
III  '  4( 1-v) 


(47) 


where  Vj  is  an  opening  displacement,  -  an  in-plane  displacement  component 
normal  to  the  crack  front,  and  Vjjj  -  an  in-plane  displacement  tangent  to  the 
crack  front.  The  above  displacements  are  for  the  nodes  closest  to  the  crack 
front,  their  distance  to  the  crack  front  being  d. 

5.4  Results  of  Computations  Using  Original  Algorithm 

The  computational  algorithm  presented  in  the  preceding  section  was  used 
to  solve  the  problem  of  the  penny  shaped  crack  under  two  kinds  of  loads: 
constant  pressure  and  constant  shear.  Three  different  meshes  were  used 
(Figure  28).  All  the  meshes  used  quarter  point  element  modeling  for  the 
displacements  close  to  the  crack  front  (which  results  in  0(/r)  variation  of 
the  displacements  as  a  function  of  distance  r  to  crack  front).  It  is 
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important  to  point  out,  though,  that  the  associated  displacement  gradient  M- 
mesh  did  not  have  the  0(7-)  singularity.  The  resulting  crack  opening 
displacement  patterns  in  the  radial  direction  along  with  the  exact  Sneddon 
[34]  and  Segedin  [36]  solutions  are  given  in  Figure  29  (z  displacement  due  to 

pressure  loading)  and  Figure  30  (x  displacement  due  to  xz  shear  loading). 

It  is  clear  from  the  figures  that  the  results  depend  very  strongly  on  the 

mesh  used  and  exhibit  a  strong  oscillatory  character.  The  accuracy  of  results 

is  clearly  unacceptable  (17?  to  25?  error  in  all  stress  concentration  factors, 
similar  order  of  magnitude  for  maximum  displacement  errors).  It  is 
interesting  to  observe,  relatively,  that  the  best  results  were  obtained  for 
the  most  uniform  mesh:  the  smoothest  of  the  three  displacement  variations  and 
17?  error  in  stress  concentration  factors. 

Stresses  in  the  crack  plane  were  calculated  for  all  three  solutions. 
Consistently,  they  showed  the  same  order  of  accuracy  as  the  crack  opening 
displacements.  Again,  the  best  results  were  obtained  for  the  uniform  mesh; 
the  errors  displayed  the  same  behavior  as  crack  opening  displacements  -  the 
biggest  errors  were  nearest  crack  front,  the  smallest  away  from  the  crack. 

One  very  clear  feature  of  all  results  was  observed:  consistently  low 
displacements  and  stresses  (as  compared  with  exact  solutions).  The  reason  for 
this  behavior  was  ascribed  to  finite  (not  zero)  compliance  of  crack  front. 

One  way  of  assuring  more  physically  appropriate  behavior  of  the  numerical 
system  was  to  implement  the  singularity  of  displacement  gradients  in  the 
vicinity  of  crack  front.  Intuitive  reasoning  would  suggest  this  as  a  way  of 
building  in  more  compliance  at  the  crack  front,  increasing  the  overall  level 
of  displacements  and  stresses. 
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Figure  29.  Crack  Opening  Displacements  Due  to  Pressure  Loading 
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5.5  Analysis  of  Method's  Equations  Using  Exact  Solution 


Poor  results  of  the  method  provoked  extensive  investigations  of  the 
sources  of  errors.  The  simplest  check  of  the  numerical  accuracy  consisted  in 
the  substitution  of  the  obtained  displacements  into  the  system  (42)  and 
calculation  of  the  right  hand  sides.  It  revealed  no  apparent  numerical 
problem. 

Further  investigation  into  the  problem  used  a  known  (Sneddon  [34]) 
solution  for  a  penny  shaped  crack  to  check  the  accuracy  of  various  matrices 
used  in  the  numerical  method.  First,  the  accuracy  of  the  interpolation  matrix 
MNM  was  checked  by  post-multiplying  it  by  exact  displacement  pattern  to  obtain 
"semi-exact"  pattern  of  displacement  gradients  (48) 


{v .  } 

1  i,aJs.e. 


=  iMNMHv1lexaot 


(48) 


Comparison  of  the  se mi-exact  displacement  gradients  with  exact  ones 
predictably  showed  the  biggest  errors,  on  the  order  of  10?,  at  the  nodes 
closest  to  the  crack  front.  Errors  at  the  very  boundary  were  obviously 
incomparable  (the  exact  solution  being  infinite  while  the  semi-exact  is 
finite)  due  to  the  simple  linear  displacement  gradient  assumption  used. 

Errors  elsewhere  on  the  crack  surface  were  small  indicating  the 

adequacy  of  the  displacement  gradient  model  for  the  regions  without 
singularity.  Overall,  the  results  of  this  test  showed  the  local  importance  of 
singularity  model  of  the  crack  front. 

The  second  test,  on  the  accuracy  of  traction  BIE  matrix  VV  (equation 
(40)),  consisted  of  substituting  noth  the  exact  and  semi-exact  displacement 
gradient  patterns  into  equation  ana  checking  resultant  right  hand 

sides.  Since  the  exact  solut.o'  . s  infinite  on  the  crack  front,  a  very  large 
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number  was  used  instead  for  the  nodes  on  the  boundary.  This  test  revealed  the 
global  importance  of  crack  front  singularity  as  all  of  the  equations  displayed 
very  large  right  hand  sides  (consistent  with  the  large  number  used  for 
infinity  at  the  crack  front). 

Ensuing  rigorous  analysis  of  the  system  of  equations  (43-43b)  proved  that 
in  fact  the  influence  of  the  crack  front  displacement  gradients  was  present  in 
all  of  the  equations  of  the  system.  Even  though  the  equations  corresponding 
to  collocation  points  on  the  crack  boundary  were  deleted  as  a  part  of 
accounting  for  boundary  conditions,  the  information  about  the  magnitude  of 
crack  front  displacement  was  present  throughout  the  system.  The  experiment 
with  the  semi-exact  displacement  gradients  showed  comparable  order  of  errors 
of  right-hand  sides  as  displacement  gradients  themselves.  The  most  important 
conclusion  of  this  test  was  the  global  importance  of  crack  front  singularity 
modeling. 

5.6  Crack  Front  Singularity  Implementation 

It  is  very  well  known  that  displacement  gradients  are  singular  at  the 
crack  front  much  as  the  stresses  are  singular  ahead  of  the  crack  (order  of 
singularity  0(^p)).  This  kind  of  singularity  may  be  modeled  effectively  by 
so-called  quarter-point  elements,  quadratic  isoparametric  finite  elements  with 
mid-side  nodes  moved  to  a  quarter  side  length  location.  In  the  current 
application,  these  elements  may  be  used  only  for  displacement  modeling  since 
the  nodal  variables  for  these  elements  are  displacements  and  the  singularity 
is  desired  in  displacement  gradients.  For  modeling  of  the  displacement 
gradients,  the  explicit  singularity  had  to  be  used  in  conjunction  with 
slightly  modified  linear  shape  functions. 

If  we  consider  a  one-dimensional  quarter  point  element,  the  displacement 
derivative  singularity  has  the  form  (Henshell  and  Shaw  [36]). 
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This  expansion  allows  for  very  accurate  modeling  of  the  displacement  gradient 
pattern  resulting  from  the  elliptical  displacement  distribution.  The  same 
variation  of  displacement  gradients  was  thus  assumed  for  both  pairs  of  4-noded 
elements.  The  resulting  shape  functions  for  these  elements  had  the  form  given 
in  eqs.  (50)  and  (51): 

Element  adjacent  to  crack  front: 


/l-1 


M2(^2}  = 


M-U^g)  =  2  (1_5r 


MJ*U152)  =  2  (1+C1) 


Element  of  the  second  row: 


M1(5152)  =  5-3{;2 


M2(«152) 


5-352  (1’el)(1‘2 


m3(c1c2) 
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In  the  above,  an  assumption  is  used  that  the  quarter  point  elements  are 
employed  for  the  first  row  of  elements  adjacent  to  crack  front.  Since  the 
shape  functions  (50)  are  unbounded  and  both  sets  are  nonlinear,  the  standard 
Galerkin  procedure  for  calculation  of  nodal  displacement  gradients  was  proven 
to  be  ineffective  and  the  collocation  method  was  used  instead  (only  on  quarter 
point  elements). 

The  effectiveness  of  the  new,  mixed  method  for  displacement  gradient 
interpolation  was  tested  on  the  mesh  of  Figure  31.  The  elliptical 
displacement  pattern  was  imposed  on  the  nodes  of  the  mesh  and  the  resulting 
displacement  gradients  were  calculated  using  standard  (non-singular)  and  new, 
mixed  method.  The  plots  of  the  results  of  both  methods,  along  with  the  exact 
and  unsmoothed  N-mesh  distributions,  are  presented  in  Figure  32.  The  typical 
behavior  of  the  interpolation  schemes  is  clearly  more  visible  in  Figure  33, 
where  it  is  shown  blown  up  on  two  elements  closest  to  the  crack  front.  The 
results  clearly  indicate  the  improved  accuracy  of  displacement  gradient 
modeling  of  the  new  mixed  method  over  the  previous,  linear  only  method.  The 
deficiency  of  the  standard  F.E.  modeling  of  displacement  gradients  is  also 
obvious  from  the  plot,  due  to  the  discontinuities  at  element  boundaries. 

Additional  improvement  in  this  displacement  gradient  modeling  is  offered 
by  the  use  of  transition  elements  (Labeyrie  and  Chauchot  [37]).  The  approach 
is  based  on  focusing  the  singularities  of  a  few  elements,  starting  from  the 
quarter  point  element,  on  a  single  location,  by  moving  their  mid-side  nodes  to 
suitable  positions  (Figure  3^).  Use  of  transition  elements  allows  for  a 
better  modeling  of  displacement  gradients  by  N-mesh  alone,  resulting  in  a 
better  approximation  by  M-mesh.  Our  study  of  the  transition  elements  showed 
that  the  better  displacement  gradient  modeling  capability  of  these  elements 
stems  from  the  presence  of  the  nonlinear  term  in  the  geometry  mapping. 


-  M-mesh  without  singularity 

-  M-mesh  with  singularity 


Use  of  Transition  Elements  in  Crack  Front 
Singularity  Modeling 


Nonuniformity  of  geometry  mapping  allows  for  nonlinear  displacement  gradients, 
enabling  a  more  accurate  fitting  of  the  crack  front  displacement  gradient 
singularity.  (Standard,  exactly  mid-side,  quadratic  elements  map  geometry  in 
a  linear  fashion.  As  a  result,  exactly  linear  variation  of  displacement 
gradients  is  possible). 

The  results  of  the  computations  for  the  previous  model,  modified  slightly 
to  produce  transition  elements,  are  presented  on  Figure  35.  They  show 
important  improvements  over  the  previous  model  (without  transition)  Figure 
32.  The  basic  difference  between  two  N-mesh  models  is  significantly  reduced 
discontinuity  of  displacement  derivatives  on  element  boundaries.  Resulting 
gains  in  the  M-mesh  interpolation  accuracy  reach  50%.  The  test  clearly 
substantiates  the  usefulness  of  transition  elements. 

The  use  of  the  new  shape  functions  for  displacement  gradient  modeling 
results  in  a  change  in  the  nodal  variables  for  the  boundary.  Instead  of  being 
displacement  derivatives  at  these  locations,  the  new  degrees  of  freedom  are 
actually  the  coefficients  A  in  the  expansion  (49): 
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where  r  is  the  crack  front  distance 

This  property  allows  for  calculation  of  stress  concentration  factors  from 
boundary  unknowns  instead  of  displacements  at  quarter  point  locations 
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5 . 7  Directions  of  Further  Work  on  Traction  BIE 

Due  to  the  lack  of  numerical  results  for  the  singularity  model  in  the 
traction  BIE,  it  is  difficult  to  formulate  definitive  conclusions.  The 
analysis  of  the  exact  solution  and  the  results  of  the  much  improved 
displacement  gradient  modeling  seem  to  indicate  the  possibility  of  substantial 
gains  in  the  stress  intensity  factor  accuracy.  A  little  less  certain  is  the 
impact  of  the  singularity  on  the  overall  level  of  displacements.  It  is 
difficult  to  estimate  the  effect  of  the  singularity  on  the  oscillatory 
character  of  displacement  results.  The  use  of  the  collocation  method  for  the 
formulation  of  system  of  equations  from  the  original  integral  equation  may  be 
responsible  for  this  phenomenon.  In  this  case,  use  of  the  weighted  residual 
approach  (Galerkin  scheme)  may  be  necessary  to  remove  the  oscillations.  The 
use  of  the  Rayleigh-Ritz  procedure  (which  is  normally  equivalent  to  the 
Galerkin  scheme)  suggested  by  Bui  (30),  in  conjunction  with  singularity 
modeling,  is  impossible,  due  to  the  appearance  of  the  (p)  term  in  the  four¬ 
fold  integral  resulting  from  the  method.  This  term  will  result  from 
two  (— )  terms  appearing  in  explicitly  formulated  singularities  of  shape 
functions . 

Further  work  in  the  development  of  the  method  will  require  studies  in  the 
method's  capacity  for  modeling  more  complex  cases.  This  capability  is 
essential  for  the  use  of  the  algorithm  as  a  part  of  alternating  method.  The 
more  complex  loading  will  probably  necessitate  the  use  of  Galerkin  scheme, 
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since  the  collocation  method  carries  no  information  about  the  interpolation  of 
loads,  into  the  final  system  of  equations.  There  are  also  indications,  Hanson 
and  Phillips  [38],  that  the  Galerkin  scheme  will  yield  a  more  tractable  set  of 
equations . 
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6.0  PROGRAM  ACCOMPLISHMENTS  AND  CONCLUSIONS 


The  principal  goal  of  the  sponsored  research  was  the  development  and 
exploitation  of  the  boundary-integral  equation  method  for  fracture  mechanics 
analysis  of  two-dimensional  problems  with  crack  tip  plasticity.  The 
formulation  approach  used  was  quite  unique  and  was  based  on  a  special  BIE 
formulation  which  explicitly  accounts  for  the  presence  of  the  crack  in 
satisfying  boundary  conditions  and,  for  the  elastic  problem,  satisfying  the 
internal  equilibrium  conditions.  Elastoplastic  behavior  was  successfully 
added  to  this  special  fracture  mechanics  modeling  capability  with  positive 
benefits  in  terms  of  numerical  and  theoretical  results. 

New  theoretical  developments  were  achieved  as  a  direct  result  of  the  use 
of  a  EIE  formulation  approach.  Finite  element  methods  begin  with  an 


approximation  of  the  total  field  equations  for  the  problem;  BIE  does  not. 
Because  the  BIE  formulation  makes  direct  use  of  relations  that  satisfy  the 
elastic  field  behavior,  even  for  the  numerical  solution  results,  certain  field 
results  are  directly  and  exactly  accounted  for: 

1.  The  analytical  formulation  established  limits  on  the  strength  of  the 
plastic  strain  increment  singularity.  In  finite  element  modeling,  one  can 
deduce  only  the  limits  on  strain  energy.  The  BIE  result  confirmed  that  the 
plastic  strain  increment  must  satisfy  the  conditions  of  singularity  previously 
adduced  to  plasticity  theory,  based  on  nonlinear  elastic  modeling. 

2.  The  analytical  result  demonstrated  that  the  plastic  singularity 
behavior  derives  from  a  singular  eigenvalue  problem  of  a  unique  form,  for  the 


plastic  field  near  the  crack  tip.  At  this  point,  the  solution  to  this 
singularity  formulation  is  imbedded  in  the  numerical  algorithm,  as  an 
analytical  solution  does  not  seem  possibj.''. 
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3.  A  new  result  was  achieved  for  the  direct  calculation  of  elastic 
stress  intensity  factors  for  cracks  influenced  by  nonsingular  inelastic 
strains.  These  problems  include  residual  strain  due  to  welding  or 
elastoplastic  notch  behavior,  thermal  strains,  etc.  The  new  result  shows  that 
extremely  accurate  stress  intensity  factor  calculations  can  be  made  which 
directly  account  for  these  volumetric  strain  distributions.  Published 
numerical  results  were  generated  for  plastic  notches  and  for  a  simplified 
welding  problem. 

4.  As  a  result  of  3.,  potentially  very  important  insight  has  been 
achieved  into  the  behavior  of  crack  extension  into  prior  plastic  zones,  due  to 
cyclic  loading  of  the  crack.  In  particular,  it  was  proven  analytically  and 
confirmed  numerically  that  the  elastic  singularity  applies  for  the  effect  of 
all  prior  crack  tip  plasticity  on  a  crack  which  has  slightly  extended  into 
this  prior  crack  tip  plasticity.  This  result  explains  why  certain  residual 
stress  retardation  models  apply  to  the  fatigue  crack  overload  problem.  In 
fact,  results  obtained  using  the  new  code,  derived  since  the  close  of  the 
current  contract,  show  that  the  effects  of  closure  and  retardation  share  a 
dual  role  on  reducing  crack  driving  force;  they  are  essentially  two  versions 
of  the. same  phenomenon. 

5.  The  last  major  analytical  result  for  the  2D  problems  was  the 
demonstration  that  the  new  BIE  code  could  be  used  as  a  direct  means  for 
calculation  of  two-dimensional  crack  weight  functions  for  problems  with 
completely  arbitrary  geometries,  boundary  conditions,  and  internal  strain 
distributions.  The  new  algorithm  was  demonstrated  for  several  test  oroblems 
and  shows  accuracy  comparable  to  the  best  numerical  models  available. 

6.  Finally,  the  results  of  the  research  will  appear  in  several  archival 
journal  articles  and  have  been  presented  in  a  variety  of  technical  symposia. 
The  foliowing  lists  the  articles  and  presentations : 
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1.  "Fracture  Mechanics,"  T.A.  Cruse,  Boundary  Element  Methods  in 
Mechanics .  book  in  series  Computational  Methods  in  Mechanics,  edited  by  D.E. 
Beskos,  Elsevier  Science  Publishers  B.V.,  Amsterdam,  to  be  published. 

2.  "Elastoplastic  BIE  Analysis  of  Cracked  Plates  and  Related  Problems, 
Part  1:  Formulation,"  T.A.  Cruse  and  E.Z.Polch,  International  Journal  for 
Numerical  Methods  in  Engineering,  vol.  23,  pp.  429-437  (1986). 

3.  "Elastoplastic  BIE  Analysis  of  Cracked  Plates  and  Related  Problems, 
Part  2:  Numerical  Results,"  T.A.  Cruse  and  E.Z.  Polch,  International  Journal 
for  Numerical  Methods  in  Engineering,  vol.  23,  pp.  439-452  (1986). 

4.  "Advanced  Algorithms  for  Fracture  Mechanics  in  Two  and  Three 
Dimensions,"  T.A.  Cruse  and  E.Z.  Polch,  2nd  International  Conference  on 
Variational  Methods  in  Engineering,  Southampton,  England,  July  17-19,  1985. 


>  w 

>  A/ 

i»aa 

y-\y. 

>Vv 
«  v  v " 


*,* 


5.  "BIE  Analysis  of  Crack  Tip  Plastic  Zones,"  T.A.  Cruse  and  E.Z.  Polch, 
Proc.  of  AIAA/ASME/ASCE/AHS  26th  Structures,  Structural  Dynamics,  and 
Materials  Conference,  Orlando,  Florida,  April  15-17,  1985,  AIAA,  New  York. 

6.  "Application  of  an  Elastoplastic  Boundary  Element  Method  to  Some 
Fracture  Mechanics  Problems,"  T.A.  Cruse  and  E.Z.  Polch,  Engineering 
Mechanics,  vol.  23,  No.  6,  pp.  1085-1096  (1986). 

7.  "A  General  Solution  Procedure  for  Fracture  Mechanics  Weight  Function 
Evaluation  Based  on  the  Boundary  Element  Method,"  T.A.  Cruse,  submitted  to 
Computational  Mechanics:  An  International  Journal. 

8.  "Buried  Crack  Analysis  with  an  Advanced  Traction  BIE  Algorithm,"  E.Z. 
Polch,  T.A.  Cruse,  and  C.-J.  Huang,  Advanced  Topics  in  Boundary  Element 
Analysis,  Ed.  by  T.A.  Cruse,  A. B .  Pifko,  and  H.  Armen,  AMD  -  Vol.  72,  (Proc. 
of  Symposium  at  ASME  Winter  Annual  Meeting,  Miami  Beach,  Florida,  November  17- 
22,  1985),  pp.  173-188,  ASME,  New  York  (1985). 

The  results  developed  to  date  are  now  being  used  in  a  systematic  study  of 
cyclic  crack  extension.  In  particular,  these  results  will  be  used  to  guide 
the  modeling  and  interpretation  of  data  for  the  current  AFOSR  sponsored 
research  into  small  crack  behavior,  Contract  F49620-84-C-0042 .  It  appears 
clear  to  us  that  the  behavior  of  fatigue  cracks  contains  much  new  information 
that  we  are  now  able  to  deduce  through  the  unique  modeling  approach  associated 
with  the  BIE  methodology. 

The  principal  focus  of  the  ongoing  work  is  the  interaction  of  prior 
plasticity  with  the  current  crack  driving  force.  Simplified  means  for 
accounting  for  this  interaction  are  needed  in  order  to  do  systematic  studies 
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of  three-dimensional  surface  cracks.  The  insight  gained  and  algorithms 
developed  in  the  first  two  year  program  form  a’  strong  basis  for  the  continued 
study  of  the  fatigue  crack  growth  problem. 
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